function u0 = twogiditer(A,F,B,G,x,epsilon,max_it,tol)
%function u0 = lagrange22(A,F,B,G)
%u0 is the result of an iterated solution with matrices of reduced size
%u0=[A;B]\[F;G];
if nargin == 4
    epsilon=1.0E-4;
    max_it = 10;
    tol = 1.0E-5;
end
[p,n] = size(B);
cvg = 0;
it_count = 0;
oldErr = inf;
newErr = 10000; 
L=1/epsilon*B'*B;
bg=1/epsilon*B'*G;
u0=(A+L)\(A*x+bg);
max_it=100;
% u1 = bicg(A+L,A*u0+bg,tol,500);

while ~cvg & (it_count <= max_it)
  u1=(A+L)\(A*u0+bg);
  %disp(['The error between iterations is ',num2str(norm(u1-u0,inf)), '.'])
  oldErr = newErr;
  newErr = norm(u1-u0,inf);%/norm(u1,inf);
  cvg = newErr <=tol;
  it_count = it_count + 1;
  u0 = u1;
end;
fprintf('Number of linear iterations is %d. ',it_count);
